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ABSTRACT 

In spherical galaxies, binary supermassive black holes (SMBHs) have difficulty reaching sub-parsec 
separations due to depletion of stars on orbits that intersect the massive binary - the "final parsec 
problem." Galaxies that form via major mergers are substantially nonspherical, and it has been 
argued that the centrophilic orbits in triaxial galaxies might provide stars to the massive binary at a 
high enough rate to avoid stalling. Here we test that idea by carrying out fully self-consistent merger 
simulations of galaxies containing central SMBHs. We find hardening rates of the massive binaries 
that are indeed much higher than in spherical models, and essentially independent of the number of 
particles used in the simulations. Binary eccentricities remain high throughout the simulations. Our 
results constitute a fully stellar-dynamical solution to the final parsec problem and imply a potentially 
high rate of events for low-frequency gravitational wave detectors like LISA. 

Subject headings: Stellar dynamics - black hole physics - 
center. 



Galaxies: kinematics and dynamics - Galaxy: 



1. INTRODUCTION 

The prospect that low-frequency gravitational radia- 
tion will be det ected by LISA , the laser-interferomete r 
space antenna (Hughes^ 120031 : iBarack fc Cutleii I2004D . 
has motivated theoretical studies into the formation and 
evolution of binary supermassive black holes (SMBHs). 
Such binaries would constitute the highest signal-to- 
noise ratio sources of low-frequency gravitational waves, 
but as is the case for virtually all potential LISA 
sources, the event rate is still poorly known, with es- 
timates r anging from a few to a few thousand events per 
year (e.g. JWvithe k Loebl[200l iRhook fc Wvithd 120051: 
ISesand 12011 " While astronomical evidence for the ex- 
istence of SMBHs with masses M, > IO^Mq in galaxy 
spheroi ds beyond the Local G roup is increasingly com- 
pelhng (|Ferrarese fc Fordll2005l ). the evidence for SMBHs 
with masses < IO^Mq, which lie more nearly in the 
LISA band, is still largely circumstantial ()van der Marell 

f)04l), and the evidence for binary SMBHs is weaker still 
(omossa 2006). 

When two galaxies merge, a binary SMBH for ms at 
the c enter of the new galaxy (|Begelman ct al. 198Qj; iRoosI 
fl98ll) . A common practice when estimating event rates 
for LISA is to equate the binary SMBH coalescence rate 
with the galaxy merger rate, the latter derived from mod- 
els of structure formation in which ga laxies me rge hierar- 
chical lv (.Haehnclt 1994; Menou, Hainian fc Naravanani 
[200lt rVolonteri. Haardt fc Madaull2003l : iWvithe' fc Loebl 



120031: iJaffe fc Backe r"2003V However it is also possible 
that coalescence is delayed, perhaps indefinitely, if the 
massive binary is limited in its ability to exchange an- 
gular momentum with surrounding stars and gas. This 
is so metimes referred to as the " final parsec problem" 
(e.g. iMilosavli evic fc Merrittll20Q3| ) since at this approx- 
imate separation, a massive binary in a spherical galaxy 
would efficiently eject a l l stars on intersect ing orbits 
(|Makino fc Funatol ' l200l iBerczik et all [2005[) . Further 
evolution of the binary is possible, even in the absence of 
gas, if these orbits are replenished, or if the galaxy geom- 
etry allows for a much larger pop ulation of centrophilic 
orbits than in the spherical case (jMerritt fc Poonll200l 
IBerczik et al|[2005l) . 

In this paper, we argue that the latter situation may 
be generic. We carry out direct A'^-body merger simu- 
lations of galaxies containing central SMBHs, and show 
that the binary SMBHs formed in such simulations con- 
tinue to shrink, via interaction with stars, well beyond 
the point at which they would be expected to stall 
in spherical galaxies. Of course, this mechanism can 
act in addition to torques from gas, if present, (e 



lEscala et al1l2005l : iDotti et al.ll2007flCuadra et al.ll200E 

or to loss-cone repopulatiq n due to "massive perturbers" 
ijPerets fc Alexanderll2008f ). etc. However it appears that 
collisionless stellar dynamics is sufficient. 

In §2 we describe the numerical methods and initial 
conditions. Simulations of isolated (spherical) models are 
presented in §3, while §4 describes the results of galaxy 
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merger simulations. §5 estimates time scales for coales- 
cence of massive binaries, and §6 sums up. 

2. NUMERICAL METHODS AND INITIAL CONDITIONS 

Our initial conditions are based on spherical galaxy 
models following Dehnen's (1993) density law: 



p{r) 



(3 - 7)Mgai ro 



and with cumulative mass profile 
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(2) 



where Mgai is the total mass of the galaxy, a is a scale 
radius, and 7 defines the inner density profile slope. We 
adopted 7 = 1 , corresponding to a Hernquist model 
(jHernquist] 1 1 9901 ) . for all initial models. A massive par- 
ticle representing a SMBH was placed at the center of 
each galaxy model. The isotropic distribution function 
reproducing p{r) in the combined potential of the stars 
and the SMBH in dynamical equilibrium was computed 
numerically and used to generate Monte-Carlo positions 
and velocities for the stars. In what follows, we adopt 
"model units" such that Mgai — G = rg = 1. 

The models can be scaled to physical units by fixing 
galaxy mass Mgai and scale radius tq combined with the 
relations 
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In a Hernquist model the effective (projected half-mass) 
radius is ~ l.Slro. For a typical, luminous (Mb ~ 
—20) elliptical galaxy or bulge, ro « 1.5 kpc, Mgai ~ 
10 "M© and [T] w 1.1 Myr. Effective radii scale very 
weakly with galaxy luminosity (though with considerable 
scatter; e.g. [Ferrarcsc et al. (2006)) and in the case of 
IO^Mq galaxy, whose SMBH would have a mass more 
suited to detection by LISA, the physical unit of time 
would be closer to 6Myr. In the latter case, the length 
of our longest integrations amounts to ~ 1.5 Gyr. 

Two sets of simulations were carried out. In the first 
set, a single spherical galaxy was created, with a central 
SMBH of mass M. = 0.001 (Model A) or 0.01 (Model B). 
(In what follows, Af, is always used to denote the mass 
of a single SMBH.) A second SMBH, of the same mass, 
was then placed at a distance of 0.5 from the center on a 
circular orbit. Integrations were continued until a time 

imax = 100. 

In the second set, two identical spherical models were 
created, each containing a single central SMBH, and the 
two galaxies were merged. Here again, two values for Af, 
were used: M. = 0.001 (Model C) or 0.01 (Model D) and 
the integrations were continued until a time imax = 250. 



TABLE 1 

Parameters of the single-galaxy integrations 



Run 


N 


A/. 


Run 


N 


M. 


Al 


32k 


0.001 


Bl 


32k 


0.01 


A2 


64k 


0.001 


B2 


64k 


0.01 


A3 


128k 


0.001 


B3 


128k 


0.01 


A4 


256k 


0.001 


B4 


256k 


0.01 


A5 


512k 


0.001 


B5 


512k 


0.01 



In order to test the dependence of the results on par- 
ticle number N , all simulations were carried out with 
TV = (32, 64, 128, 256, 512) x 10^. Tables [Hand H summa- 
rize the parameters of the two sets of runs. 

The A^-body integrations we re carried out us- 
ing (?!)GRAPE (|Harfst et al.l I2007D . a parallel, direct- 
summation A^-body code that uses special-purpose hard- 
ware to compute pairwise forces between all particles: 



N 



mj{vi-Y,) 



2^3/2 ■ 



(4) 



Here nii and ri are the mass and position of the ith 
particle respectively and e is a force softening param- 
eter. We set e = 10~^ in all the runs. ^GRAPE in- 
tegrates the equations of motion using a fourth-order 
Hermite integrator with individual block time steps 
simila r to the well established Aarseth direct N-body 
codes ljMakino fc Aarsethlll992| ). The time step of parti- 
cle i at time t was computed from the standard formula: 



At; = 



|a(0||a^(t)| + |aH^)P 
|ai(t)||a3(t)| + |a2(t)| 



(5) 



where a is the acceleration of ith particle and the super- 
script fc indicates the order of the time derivative. The 
time-step parameter r\ was set to 0.01 for all our runs. 
The ch oice of parameters was motivated bv lBerczik et al.) 
()2005[) who found that integration accuracy is more sen- 
sitive to T] than to e. Those authors adopted 77 = 0.01 
and e = 1 X 10^^ for their study of SMBH binary evolu- 
tion. In order to resove the orbit of the massive binary 
in our simulations down to a semi-major axis a « 10"** 
we adopted the smaller softening length e = 10~^. With 
these choices of 77 and e, the relative energy error in our 
simulations was ~ 10^^. For a more detailed discus- 
sion of the se issues we refer the reader to section 2 of 
IMerritt et a l. (2007) . 

The A-body integrations were carried out on three, 
special-purpose computer clusters. Titan, at the 
Astronomisches Rechen-Institut in Heidelberg, and 
gravitySimulator, at the Rochester Institute of Tech- 
nology are both 3 2-node clusters em ploying GRAPE 
accelerator boards ijHarfst et al.ll2007t ). Some calcula- 
tions were also carried out on the GPU-enhanced cluster 
"Kolob" at the University of Heidelberg. For the latter 
runs we used the modified versio n of the (/iGRAPE cod e 
including the SAPPORO library ()Gaburov et al.|[2009l) . 

3. ISOLATED MODELS 

Evolution of a binary SMB H in a gas-free galaxy can 
be divided into three phases (|Begelman et aliri980[ ). In 
first phase, the SMBHs are unbound to each other and 
move independently in the galaxy potential. The evolu- 
tion of the individual SMBH orbits during this phase is 
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Fig. 1. — Evolution of the separation between the SMBHs, in 
ten A'-body integrations of isolated spherical galaxies with different 
particle numbers according to Table [T] Models A (top) have M, = 
0.001 and models B (bottom) have M, = 0.01. 
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Fig. 2. — Evolution of the inverse semi-major axis of the massive 
binary in the isolated models. 



governed by dynamical friction against the background 
stars. The second phase begins when the two SMBHs 
are close enough together to form a bound pair. In this 
phase, the separation between the SMBHs decreases very 
rapidly due to the combined effects of dynamical friction 
and ejection of stars by the gravitational slingshot. (For 
a detailed discussion of this phase, we refer the reader to 
the careful analysis in iMilosavlievic fc Merrit"^ ()2001|) .) 
At the end of the second phase, the two SMBHs form a 
hard binary, defined as a binary that ejects passing stars 
with positive (unbound) energies. Phase three consists of 
slow evolution of the binary as depleted orbits are repop- 
ulated, via two-body scattering or some other process. 
In a spherical galaxy with large N, the timescale for or- 
bital re-population is essentially the two-body relaxation 
time which scales as ^ N/ \ nN. This is the origin of t he 
"final-parsec problem" ((Merritt k Milosavlievicll2005D . 

Figure [T] plots the separation R between the two SMIBH 
particles versus time in the isolated galaxies. Models A 
and B. Here we can see the three phases of binary evo- 
lution described above. 

(1) In the first phase, the separation between the two 
SMBHs decreases due to dynamical friction. This phase 
ends when i? r^, the gravitational influence radius, 
can be defined as the radius of a sphere around the binary 
SMBH that encloses a stellar mass equal to twice the 
binary mass. For Models A, M, = 0.001 and rh « 0.07. 



From FigurelU we can see that the first phase ends indeed 
when R « 0.07. Here we want to emphasize that this 
estimate does not take into account the changes in the 
galaxy that are induced by the presence of the second 
SMBH. 

For Models B, the radius of influence is ^ 0.25. Judg- 
ing from the bottom panel of Figure[Tl it appears that the 
flrst phase ends approximately when R w 0.25, although 
the transition is not as sharp as in model A, presumably 
because the initial separation is only ~ 2r/i. 

(2) In the second phase, the separation between the 
two SMBHs decreases very rapidly. Dynamical friction, 
and the ejection of stars by the gravitational slingshot, 
act together to efficiently extract angular momentum 
from the massive binary. These processes are not well 
defined in the regime where the massive binary is nei- 
ther very hard nor v ery soft; equations (13) and (14) of 
IMilosavlievic fc Mer ritt (2001) give approximate expres- 
sions for the rate at which energy is transferred to stars 
by the two mechanisms. For Models A this process con- 
tinues until t « 40 while for Models B, this phase ends 
at i w 20. The motion of two SMBHs in this phase is 
approximately Keplerian. We define the semi-major axis 
a and eccentricity e via the standard relations: 
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a 
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(6) 
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e = Wl 



2h^ 
1^ 



R 



(7) 



where ^ — G(A/,i + Af,2), v is the relative speed, and h 
is the specific angular momentum of the relative motion. 

(3) The rapid phase of binary hardening comes to an 
end when a Ri ah, where ah is the semi-major axis of a 
"hard binary" : 

(e.g. iMerritt et all '2007). Here q = M.2/Af.i) is the 

mass ratio of two SMBHs; for both Models A and B, 
q=l. For Models A, ah ~ 0.004 with a^^ w 250 and for 
Models B, ah « 0.016 with « 65. 

From Figure [2l we see that there is no clear depen- 
dence of the binary hardening rate on N until the bi- 
nary becomes hard, a < ah- For Models A this hap- 
pens around t = 25 and for Models B around t = 
30. Beyond these times the binary hardening rate ex- 
hibits a clear N dependence in the sense that harden- 
ing is slower for larger N. Similar A^-dependence has 
been seen in othe r studies of binary evolution in spher- 
ical galaxies (e.g. iMakino &: Funatol 120041 : iBerczik et al.l 
I2006t IMerritt et al.l 120071 ). For real spherical galaxies, 
with much larger N, the binary would stop evolving be- 
yond this point. 

We estimated the hardening rate s in the A"-body in- 
tegrations, where 



(9) 



dt \ a 



by fitting straight lines to a ^ (t) in an interval At = 50 
from t = 50 to t = 100. Figure [3] shows the N- 
dependence of s. We see that for A^ > 50k, the hardening 
rate is a decreasing function of A^; scaling as oc A^~°-^; 
this is consistent with the scalin gs found by other au - 
thors in the same A'^-range, e.g. iBerczik et all (pOOl : 



IMerritt etlll (120071 ). Comparison of Figures [3ji and 
b also shows lower hardening rates for Model B com- 
pared with Model A, by roughly an order of magni- 
tude. For large A^, the two-body relaxation time is long 
compared with orbital periods of stars near the massive 
binary, and the loss cone around the binary is nearly 
empty. In this regime, and in the asymptotic (large-A^) 
limit, the rate of binary evolution is predicted to scale as 
~ \(M ,i + M.a)]"^ (e.g. equation (32) of IMerritt et al.l 
((200l ). Given that a is roughly a factor ten smaller, at 
a given time, in Models A as compared with Models B, 
there is an additional delay in the early phase of the hard 
binary evolution. However this delay is negligible for the 
decay time Iq until which the semimajor axis reaches ao, 
when energy loss by GWs start to dominate. In section 
|5| a more detailed discussion of the total time to merge 
the binary SMBH is given. 

Figure [5] shows the density profile for model A5 in the 
final phase. The central logarithmic slope of the density 
profile has dropped from —1 to ^ —0.6. 

4. GALAXY MERGERS 

We studied the evolution of binary SMBHs in merging 
galaxies by creating two identical galaxy models, each 
with a central SMBH, and placing the two galaxies on 
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Fig. 3. — Af-dependence of the binary hardening rate in the 
isolated galaxy models. Top: Models A; bottom: Models B. 



TABLE 2 

Parameters of the galaxy merger simulations 



Run 


2 X Af 


M. 


CI 


64k 


0.001 


C2 


128k 


0.001 


C3 


256k 


0.001 


C4 


512k 


0.001 


Dl 


64k 


0.01 


D2 


128k 


0.01 


D3 


256k 


0.01 


D4 


512k 


0.01 



bound relative orbits. Table [5] summarizes the parame- 
ters of the galaxy merger models. The initial separation 
of the two galaxy centers was 20ro. The initial relative 
velocity of the two galaxies was chosen such that the 
SMBH separation at first pericenter passage was ~ rp; 
in other words, the initial orbit of the binary galaxy was 
substantially eccentric. This was done in order to reduce 
the computational time required to bring the two SMBHs 
close together; in fact, a galaxy merger from such nearly 
"head-on" initial conditions is unlikely. 

Figure 0] illustrates the merging of two galaxies. In 
the initial phase (top panels) the two SMBH particles 
remain strongly associated with their respective density 
cusps. The two cusps merge into one a,t t ^ 100, follow- 
ing which the cusp is gradually destroyed by the gravi- 
tational slingshot ejection of stars (Fig. [5]). 
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Fig. 4. — Density contours projected onto the initial orbital 
plane for Model C3. First row: t = (0,70,90). Second row; t = 
(110,150,200). 

Figure |6] shows the separation between the two SMBHs 
in the merger models. The first pericenter passage occurs 
around t « 80. As the galaxies merge, the two SMBHs 
remain centrally located in their respective density cusps, 
until coming close enough together that a binary SMBH 
forms, at t w 100, or slightly earlier in the case of Model 
D (phase one). The separation between the two SMBHs 
then decreases very rapidly due to the combined effects of 
dynamical friction and slingshot ejection of stars (phase 
two). Once the hard binary is formed, the hardening 
rate of the binary decreases (phase three) . This behavior 
is qualitatively similar to what was seen in the isolated 
galaxy models. 

However there is one important difference between the 
binary evolution in the isolated and merging galaxies. 
Figure [7] shows the evolution of the binary semi-major 
axis, and Figure [S] shows the binary hardening rates in 
the galaxy merger simulations; the latter were computed 
in the same way as in the isolated galaxy models, by 
linear fit of 1/a in the range 150 < t < 200. In the galaxy 
merger models, there is essentially no dependence of the 
binary hardening rate on particle number. Furthermore, 
in Models C and D, the hardening rates are much higher 
- more than five times higher than in the Models A and 
B, respectively, when N = 512k. 

A likely explanation for this difference is that the non- 
spherical shapes of the merger remnants provide an ad- 
ditional source of torque on the stellar orbits, allowing 
them to interact with the central binary on a timescale 
shorter than the two-body relaxation time. It has been 
argued (jMerritt fc PoonI 120041: [Berczik et aLll2006D that 
allocating even a small fraction of the stellar orbits to 
a "centrophilic" family could lead to binary hardening 
rates that are much larger than those that result from 
two-body relaxation alone. 

We evaluated the shapes of the merged galaxies in a 
number of ways. The central density contours for Model 
C3 are shown in Figure |9l Departures from axisymmetry 
are evident. Figure [TOl shows the principle axis ratios of 
the merger remnant in Model C3 as a function of time; 
these were defined as the axis ratios of a homogeneous 
ellipsoid with the same inertia tensor. The departures 
from spherical symmetry are modest, but definite, and 
they appear to be nearly independent of time toward the 



end of the simulation. 

The evolution of binary eccentricity in the merger mod- 
els is presented in Figure [11] We notice high eccen- 
tricities for both models as soon as the SMBHs become 
bound, approaching unity for models D. The high eccen- 
tricities are due in part to the high eccentricity (e = 0.95) 
of the galaxies' relative orbit prior to the merger, and, in 
some runs, to post-merger evolution. 

5. TIME SCALES FOR COALESCENCE 

SMBH binaries are a potentially important source of 
gravitational wave (GW) emission. Gravitational waves 
extract energy and angular moment um from the b inary, 
hence changing its orbital elements. iPetersI (|1964i) gives 
approximate, orbit-averaged expressions for the rates of 
change of a binary's semimajor axis and eccentricity due 
to GW emission: 



GW 



GW 
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We define tGw(oo, cq) as the time required, according to 
the coupled equations (llOal llObp . for the binary semi- 
major axis to shrink to zero under the infiuence of GW 
emission, starting from the initial values (ao,eo). This 
time is a strong function of bq for cq ~ 1, as it is in some 
of our simulations. 

The full time to coalescence, tcoai, includes also the 
time from the start of the simulation until the GW regime 
is entered. We estimated icoai as follows. In the iV-body 
merger models, the post-merger hardening rate of the bi- 
nary, s, is essentially independent of a and N (Figures [71 
[8]). The rate of change of the binary semi-major axis in 
this regime is approximately 



da\ 



const. 



(11) 



Given values for Afgai and rp, this dimensionless rate 
can be converted into physical units. At some time 
in the simulation, a will have fallen to such a small 
value that {da/dt)-^<i^ as given by equation ()11|) will equal 
{da/dt)Qyf, where the latter quantity is computed from 
both a and e. We define to as the time at which these 
two rates are equal. The full time to coalescence, icoai, 
is the time from the start of the simulation until to i plus 
iGw(ao, eo), where oq = a{t = to), eo = e(t to)- 

In most cases, to exceeded the final time tfinai of the 
simulation. In these cases, the time At between tgnai and 
to was assumed to be 



At = S(tfi„al)^^ 



1 

ao 



1 



(12) 



In such cases, we also assumed that e(t) was constant, 
e(t) = e(tfinai), for t > tfinai whcu computing tew- 

Computation of the GW evolution rates required adop- 
tion of particular values for the physical units of length 
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Fig. 5. — Left: Spatial density profile for model A5 averaged over five different time steps (t = 60, 68, 76, 84, 92). The figure also shows 
the initial density profile (7 = 1) and a fit to the data using Eq. [T] with 7 = 0.6 and ro = 0.8. The central density drops as the binary 
ejects mass from the core. Right: Spatial density profile for model C4 averaged over five different time steps (i = 140, 160, 180, 200, 220). 
The figure also shows the initial density profile (7 = 1) and a fit to the data using equation [T] with 7 = 0.3 and ro = 0.6. 
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Fig. 6. — Evolution of the separation between the two SMBHs 
in the galaxy merger simulations. 
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Fig. 7. — Evolution of the inverse semi-major axis of the binary 
SMBH during the galaxy merger simulations. 



and time, hence of Mgai and tq. Table 3 gives resuhs 
for two different assumed values of the galaxy mass: 

Mgal = (109,10")MQ. 

Because of the strong dependence of tew on eccentric- 
ity, Table 4 gives oqw and At, for run C3 only, under 
two different assumptions about cq: = 0.5 and eo = 0. 
The total decay time increases by a factor of 2-3 if the 
orbit is assumed circular at the start of the GW phase. 



For low mass-galaxies (lO^Mo), the total time to 
merge the SMBHs exceeds ^ 1 Gyr, whereas for high- 
mass galaxies (lO^^M©) the coalescence time is well be- 
low 0.5 Gyr. This is contrary to the usual statement that 
binary coalescensce is slower in more massive binaries. 
That statement reflects the much longer timescales for 
gravitational encounters in massive galaxies; but in our 
merger models, evolution prior to the GW regime occurs 
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Fig. 8. — Binary hardening rates in the galaxy merger simula- 
tions (top panel: Model C, bottom panel: Model D). 
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Fig. 10. — Ratio of intermediate to major (b/a) and minor to 
major (c/a) axes for Run C3. 
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Fig. 9. — Density contours for merger Model C3 at t = 150, 
projected on xy plane(left), xz plane(center), yz plane(right). 



much more rapidly than predicted based on collisional 
loss cone repopulation rates, and the dominant influence 
on the coalescence time is the timescale for GW energy 
loss, which is longer for smaller SMBH masses. We note 
that all of the coalescence times in Table 3 are short 
compared with those expected in spherical models (e.g. 
Figure 19 of Merritt, Mikkola & SzeU 2007). 

The lower mass scaling corresponds to sources de- 
tectable by LISA. We find coalescence times of a few Gyr 
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Fig. 11. — Evolution of binary eccentricity for models C (top) 
and models D (bottom). 

in these cases, which, while relatively short, are never- 
theless long enough that SMBH binaries in these galax- 
ies may not achieve full coalescence before a subsequent 
galaxy merger occurs. 

6. SUMMARY 

The "final parsec problem" refers to the difficulty of 
bringing binary supermassive black holes (SMBHs) to 
full coalescence, due to rapid depletion, via slingshot 
ejection, of stars on orbits that intersect the massive 
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binary. It has been suggested ([Merritt fc PoonI 120041) 
that binary hardening can occur much more efficiently 
in triaxial galaxies due to the greater population of cen- 
trophilic orbits. Generically, galaxies formed in dissi- 
pationless mergers are triaxial in shape, and one might 
expect to avoid the final parsec problem in the case of 
binaries that form as a result of galaxy mergers. To test 
this hypothesis, we carried out high-accuracy, direct TV- 
body simulations of dissipationless mergers containing 
SMBHs. The results were compared with a second set of 
simulations of massive binaries in spherical models. The 
strong A'^-dependence of the binary hardening rate in the 
spherical geometry was found to vanish in the merger 
models. This was true, even though the departures from 
axisymmetry in the merged galaxies were slight. The bi- 
nary SMBHs in our merger simulations formed with large 
eccentricities, due to the merger geometry, and eccen- 
tricities remained high during the subsequent evolution 
of the binary. Our results suggest that the "final parsec 
problem" may not be as serious as previously believed, 
and that prompt coalescence of binary SMBHs following 
mergers may be common, even in galaxies lacking gas. 



This work began as a student project during the Ad- 
vanced School and Workshop on Computational Gravita- 
tional Dynamics held at the Lorentz Center, Leiden from 
3 May 2010 - 13 May 2010. We thank Simon Portegies 
Zwart for his efforts in organizing this workshop. We 
thank Peter Berczik and Ingo Berentzen for their sup- 
port on the numerical simulations and helpful discus- 
sions. FK was supported by a grant from the Higher Ed- 
ucation Commission (HEC) of Pakistan administrated by 
the Deutscher Akademischer Austauschdienst (DAAD). 
DM was supported by grants AST-0807910 (NSF) and 
NNX07AH15G (NASA). This research and the computer 
hardware used in Heidelberg were supported by project 
"GRACE" 1/80 041-043 of the Volkswagen Foundation, 
by the Ministry of Science, Research and the Arts of 
Baden- Wlirttcmberg (Az: 823.219-439/30 and 823.219- 
439/36), and in part by the German Science Founda- 
tion (DFG) under SFB 439 (sub-project Bll) "Galaxies 
in the Young Universe" at the University of Heidelberg. 
The HPC - GPU cluster "Kolob" is supported by the 
University of Heidelberg through the excellence initiative 
for German universities by the innovation funds FRON- 
TIER. 



REFERENCES 



Barack, L., & Cutler, C. 2004, Physical Review D, 69, 082005 
Begelman, M. C, Blandford, R. D., & Rees, M. J. 1980, Nature, 
287, 307 

Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680 
Berczik, P., Merritt, D., Spurzem, R., Bischof, H., 2006, ApJ, 642, 
L21 

Cuadra, J., Armitage, P. J., Alexander, R. D., & Begelman, 

M. C. 2009, MNRAS, 393, 1423 
Dehnen W., 1993, MNRAS, 265, 250 

Dotti, M., Colpi, M., Haardt, P., & Mayer, L. 2007, MNRAS, 
379, 956 

Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2005, 

ApJ, 630, 152 
Ferrarese, L., et al. 2006, ApJS, 164, 334 

Ferrarese, L., & Ford, H. 2005, Space Science Reviews, 116, 523 
Gaburov E., Harfst S., Portegies Zwart S., 2009, New Astron., 14, 
630 

Gebhardt, K., et al. 1996, Astron. J. 112, 105 

Haehnelt, M. G. 1994, Mon. Not. R. Astron. Soc. 269, 199 

Harfst, S., Gualandris, M., Merrit D., Spurzem, Portegies Zwart 

S., Berczick P., 2007, NewA, 12, 357 
Hashimoto Y., Funato Y., Makino J., 2003, ApJ, 582, 196 
Hernquist, L., 1990, ApJ, 356, 359 

Hughes, S. A. 2002, Mon. Not. R. Astron. Soc. 331, 805 

Hughes, S. A. 2003, Annals of Physics, 303, 142 

Jaffe, A. H. & Backer, D. C. 2003, Astrophys. J. 583, 616 



Komossa, S. 2006, Memorie della Societ Astronomica Italiana, 77, 
733 

Makino, J., & Funato, Y. 2004, ApJ, 602, 93 

Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 

Makino, J., & Funato, Y. 2004, ApJ, 602, 93 

Menou, K., Haiman, Z., & Narayanan V. K. 2001, Astrophys. J. 

558, 535 
Merritt, D. 2006, ApJ, 648, 976 

Merritt, D., & Milosavljevic, M. 2005, Living Reviews in 

Relativity, 8, 8 
Merritt, D., Mikkola, S., Szell, A. 2007, ApJ, 671, 53 
Merritt, D., & Poon, M. Y. 2004, ApJ, 606, 788 
Milosavljevic, M., & Merritt, D. 2001, ApJ, 563, 34 
Milosavljevic, M., & Merritt, D. 2003, The Astrophysics of 

Gravitational Wave Sources, 686, 201 
Nakano, T., Makino, J. 1999, ApJ, 510, 155 
Nakano, T., Makino, J. 1999, ApJ, 525, L80 
Perets, H. B., & Alexander, T. 2008, ApJ, 677, 146 
Peters, P. C. 1964, Phys. Rev. B, 136, 1224 
Rhook, K. J., & Wyithe, J. S. B. 2005, MNRAS, 361, 1145 
Roos, N. 1981, Astron. Astrophys. 104, 218 
Sesana, A. 2010, ApJ, 719, 851 

van der Marel, R. P. 2004, Coevolution of Black Holes and 
Galaxies, 37 

Volonteri, M., Haardt, F. & Madau, P. 2003, Astrophys. J. 582, 
559 

Wyithe, J. S. B. & Loeb, A. 2003, Astrophys. J. 590, 691 



Efficient Merger of Binary Supermassive Black Holes in Merging Galaxies 



9 



TABLE 3 

Time to Gravitational Wave Coalescence 



Run 


"final 




final 


eo 


ao 


(pc) 




to (Gyr) 


tcoai (Gyr) 


CI 


1.2 X 10- 


-4 


120.3 


0.51 


(2.6 X 10" 


-3 


1.2 X 10" 




(11.5,0.54) 


(15.9,0.94) 


C2 


7.2 X 10" 


-5 


115.2 


0.96 


(1.2 X 10" 


-2 


8.5 X 10" 




(2.6,0.21) 


(5.1,0.29) 


C3 


7.3 X 10" 


-5 


108.3 


0.80 


(4.9 X 10" 


-3 


3.3 X 10" 




(7.5,0.35) 


(10.8,0.41) 


C4 


6.8 X 10" 


-5 


118.2 


0.95 


(1.0 X 10" 


-2 


8.3 X 10" 




(2.8,0.22) 


(5.7,0..30) 


Dl 


5.2 X 10- 


-4 


13.0 


0.97 


(9.0 X 10- 


-2 


6.7 X 10" 




(2.7,0.21) 


(5.5,0.35) 


D2 


4.0 X 10" 


-4 


11.2 


0.88 


(4.4 X 10- 


-2 


3.0 X 10- 




(6.7, 0.26) 


(10.8,0.47) 


D3 


5.1 X 10" 


-4 


10.7 


0.95 


(7.2 X 10- 


-2 


5.4 X 10- 




(4.0,0.22) 


(7.6, 0.38) 


D4 


4.7 X 10- 


-4 


11.0 


0.98 


(9.0 X 10- 


-2 


9.1 X 10- 




(2.2,0.23) 


(4.7, 0.26) 



TABLE 4 

Time to Gravitational Wave Coalescence (Run C3) 



eo 


ao (pc) 


tcoai (Gyr) 


0.5 
0.0 


(2.6 X 10-^,1.8 X 10"^) 
(1.9 X 10-3,1.3 X 10-2) 


(14.96,0.65) 
(20.46,0.86) 



